On the universality of anomalous one-dimensional heat conductivity 
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In one and two dimensions, transport coefficients may diverge in the thermodynamic iimit due to 
iong-time correlation of the corresponding currents. The effective asymptotic behaviour is addressed 
with reference to the problem of heat transport in Id crystais, modeled by chains of ciassical nonlinear 
oscillators. Extensive accurate equilibrium and nonequilibrium numerical simulations confirm that 
the finite-size thermal conductivity diverges with the system size L as k oc L a . However, the 
exponent a deviates systematically from the theoretical prediction a — 1/3 proposed in a recent 
paper [O. Narayan, S. Ramaswamy, Phys. Rev. Lett. 89, 200601 (2002)]. 
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Strong spatial constraints can significantly alter trans- 
port properties. The ultimate reason is that the re- 
sponse to external forces depends on statistical fluctu- 
ations which, in turn, crucially depend on the system 
dimensionality d. A relevant example is the anomalous 
behaviour of heat conductivity for d < 2. After the pub- 
lication of the first convincing numerical evidence of a 
diverging thermal conductivity in anharmonic chains , 
this issue attracted a renovated interest within the theo- 
retical community. A fairly complete overview is given in 
Ref. , where the effects of lattice dimensionality on the 
breakdown of Fourier's law are discussed as well. Anoma- 
lous behaviour means both a nonintegrable algebraic de- 
cay of equilibrium correlations of the heat current J(t) 
(the Green-Kubo integrand) at large times t -> oo and 
a divergence of the finite-size conductivity k(L) in the 
L — > oo limit. This is very much reminiscent of the prob- 
lem of long-time tails in fluids where, in low spatial di- 
mension, transport coefficients may not exist at all, thus 
implying a breakdown of the phenomenological constitu- 
tive laws of hydrodynamics. In Id one finds 

k{L) oc L a , (J(t)J(O)) oc , (1) 

where a > 0, — I < 5 < 0, and ( ) is the equilibrium av- 
erage. For small applied gradients, linear-response the- 
ory allows establishing a connection between the two ex- 
ponents. By assuming that k(L) can be estimated by 
cutting-off the integral in the Green-Kubo formula at the 
"transit time" L/v (v being some propagation velocity of 
excitations), one obtains k oc L~ s i.e. a = —S. 

Determining the asymptotic dependence of heat con- 
ductivity is not only important for assessing the univer- 
sality of this phenomenon, but may be also relevant for 
predicting transport properties of real materials. For in- 
stance, recent molecular dynamics results obtained with 
phenomenological carbon potentials indicate an unusu- 
ally high conductivity of single- walled nanotubes Q: a 
power-law divergence with the tube length has been ob- 



served with an exponent very close to the one obtained 
in simple Id models 

The analysis of several models Q clarified that anoma- 
lous conductivity should occur generically whenever mo- 
mentum is conserved. For lattice models, this amounts 
to requiring that at least one acoustic phonon branch is 
present in the harmonic limit. The only known exception 
is the coupled-rotor model, where normal transport is 
believed to arise as a consequence of the boundedness of 
the potential. 

It is thus natural to argue about the universality of 
the exponent a. On the one hand, there exist two theo- 
retical predictions, namely a — 2/5, which follows from 
self-consistent mode-coupling theory plla], and a = 1/3, 
obtained by Narayan and Ramaswamy [9( by performing 
a renormalization group calculation on the stochastic hy- 
drodynamic equations for a Id fluid. On the other hand, 
the available numerical data for a range from 0.25 to 
0.44. The most convincing confirmation of the 1/3- value 
have been obtained by simulating a one-dimensional gas 
of hard-point particles with alternating masses [ToL 111) 
and random-collision models In the former case, a 
careful determination of the scaling exponent is however 
hindered by the presence of large finite-size corrections 
that are still sizeable for O(10 4 ) particles. As a matter 
of fact, other authors 0] report significantly smaller es- 
timates (a ~ 0.25). This anomaly is possibly due to the 
lack of microscopic chaos in that model ^lj ■ The results 
obtained for models of Id crystals are more controversial, 
but consistently larger than 1/3. For instance, in the case 
of the Fermi-Pasta-Ulam (FPU) chain, the best estimate 
sofar is a ~ 0.37 HQ. 

However, with the only exception of Ref. 0], all nu- 
merical investigations limit themselves to fitting the scal- 
ing behaviour in a suitable range, without determining 
possible finite-size corrections, so that none of them can 
be fully trusted. In view of the general relevance of estab- 
lishing the existence of one or more universality classes, 
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in the present paper we present far more accurate sim- 
ulations which allow determining the effective exponents 
for different lengths and frequencies. We anticipate that 
a is definitely larger than 1/3 in a Id crystal model and 
possibly in agreement with the mode-coupling prediction. 

We consider an array of N point-like identical atoms 
ordered along a line. The position of the n-th atom is 
denoted with x n , while its mass is fixed, without loss 
of generality, equal to unity. By further assuming that 
interactions are restricted to nearest-neighbour pairs, the 
equations of motion write 



-F n + F n -\ 



F n 



-V'(x nA 



(2) 



where V'(z) is a shorthand notation for the first deriva- 
tive of the the interparticle potential V with respect to 
z. The microscopic expression of the heat current is 



J = E 



2^+1 



(3) 



where h n is a suitably defined local energy 0] • For small 
oscillations (compared to the lattice spacing b = L/N), 
the second term can be neglected and x„ — x n ~\ ~ b, so 
that Eq. can be approximated by 



J 



■£(* 



n+l 



(4) 



The customary way to evaluate the thermal conduc- 
tivity k is through the Green-Kubo formula 
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K-GK 



kuT 2 t—>oc 



dr lim L- x (J(r)J(0)) 



L — ►oo 



(5) 



A crucial, sometimes overlooked |l5j . point is that such 
formulae are formally identical for different statistical en- 
sembles, but the definition of J differs, because of "sys- 
tematic" contributions associated with other conserva- 
tion laws that must be subtracted out |l6l| . For instance, 
expression |@J is correct in the microcanonical ensemble 
with zero total momentum, while in the canonical ensem- 
ble (for large N) it is 

J = t| £(£n+l + x n ) Fn ~ bv (^F n ) , (6) 

n n 

vo being the center-of-mass velocity. This choice ensures 
that the autocorrelation of J vanishes for t — ► oo. 

With reference to Eq. 0), the possibly anomalous be- 
havior can be analyzed by computing the power spec- 
trum of the heat current J. Since we are interested in 
the long- wavelength and small-frequency behavior, it is 
convenient to consider a highly nonlinear model in the 
hope that the asymptotic regime sets in over shorter time 
and space scales. Moreover, it is advisable to work with a 
computationally simple expression of the force. The best 



compromise we have found is the quartic Fermi-Pasta- 
Ulam potential 



V{z) = \ {z-af 



(7) 



Indeed, after the change of coordinates x n = u n + na |l7|. 
the physical distance a disappears from the equations of 
motion for u n . This model has no free parameters: since 
the potential expression is homogeneous, the dynamics is 
invariant under coordinate rescaling, so that the energy 
per particle e can be set, without loss of generality, equal 
to 1. 

First, we have performed equilibrium microcanonical 
simulations by integrating Eqs. (J2J (with periodic bound- 
ary conditions) with a fourth-order symplectic algorithm 
jl8| . The power spectra S(f) of J are reported in Fig.^ 
The lowest curves are data for the quartic FPU model 
obtained by averaging over 30,000 random initial condi- 
tions. In order to further decrease statistical fluctuations, 
a binning of the data over contiguous frequency intervals 
has been performed as well. 




FIG. 1: Power spectra of the flux J as defined in Eq. Q. The 
lowermost curves refer to the purely quartic FPU model £J 
with N = 2048 (solid) and 1024 (dashed). The upper curves 
correspond to the repulsive FPU model - Eq. 1101 - for a = 2, 
2.3, and 2.5 (from top to bottom) and N = 1024 (solid) or 
512 (dashed). All microcanonical simulations are performed 
for the same energy density e = 1 with time step h — 0.05 for 
10 6 — 10 7 steps. For clarity, the curves have been arbitrarily 
shifted along the vertical axis. 

The long-time tail Q manifests itself as a power-law 
divergence f s in the low-/ region. By comparing the 
results obtained for different numbers of particles, one 
can clearly see that finite-size corrections are negligible 
above a size-dependent frequency f c (N). By fitting the 
data in the scaling range [f c (N), f s ], where f s ~ 10~ 3 , 
we find 5 — —0.39(6). These values are consistent with 
previous, less-accurate, findings for similar models, such 
as the standard FPU @, Q and the diatomic Toda 
chains, thus confirming the expectation that they all be- 
long to the same universality class. 

In order to perform a more stringent test of the scaling 
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FIG. 2: The logarithmic derivative 5 e fi of the energy-flux 
spectrum versus the frequency / for the pure FPU quartic 
potential - FPU-4 with N = 2048 - and for model iJTDJ with 
a = 2, a — 2.3, and a = 2.5 and N = 1024, respectively. The 
two horizontal lines correspond to the theoretical predictions 
— 1/3 and —2/5. The statistical error is on the order of the 
observed irregular fluctuations. 



behavior, we have determined the logarithmic derivative 



M/) = 



(8) 



for different frequencies. Since finite-size effects are re- 
sponsible for the saturation of S(f) when / — > 0, f c (N) 
can be identified (see Fig-El as the frequency below which 
5 e ff starts growing towards zero. Above / c , the quality 
of our numerical data allows revealing a slow but sys- 
tematic decrease of S g upon decreasing /, which ap- 
proaches —0.44, a value that is incompatible not only 
with the renormalization-group prediction of Ref.p , but 
also with the result of mode-coupling theory 0, H ■ Fur- 
thermore, convergence seems not fully achieved in the 
accessible frequency range. 

Accordingly, it is advisable to look at thermal conduc- 
tivity by means of nonequilibrium simulations too. It is 
sufficient to measure the heat flux in a system put in 
contact with two heat reservoirs operating at different 
temperatures T+ and T_. Several methods have been 
proposed based on both deterministic and stochastic al- 
gorithms 0- Regardless of the actual thermostatting 
scheme, after a transient, an off-equilibrium stationary 
state sets in, with a net heat current flowing through the 
lattice. The finite-size thermal conductivity k(L) is then 
estimated as the ratio between the average flux J and the 
overall temperature gradient (T+ — T-)/L. Notice that, 
by this definition, k amounts to an effective transport 
coefficient including both boundary and bulk scattering 
mechanisms. 

We have used the Nose-Hoover thermostats described 
in detail in Ref. 0. In order to fasten the convergence 
towards the stationary state, the initial conditions have 
been generated by thermostatting each particle to yield a 
linear temperature profile (see This method is very 
efficient, especially in long chains, when bulk thermaliza- 



tion may be significantly slow. The heat flux J has been 
obtained by integrating the equations over more than 10 6 
time units and by further averaging over a set of about 30 
initial conditions. Simulations of the quartic FPU model 
with chains of length up to 65536 sites and free bound- 
ary conditions exhibit again a systematic increase of the 
effective exponent 



a e s(L) = 



d In k 



(9) 



as it can be seen from the full dots in Fig. although 
one can also observe that the four rightmost values are in 
very good agreement with the mode-coupling exponent. 

In order to compare more closely equilibrium and 
nonequilibrium simulations, one can assume, following 
the argument exposed below Eq. Q), that the finite-size 
conductivity k(L) is determined by correlations up to 
time t = L/v s , where v s is the sound velocity. This 
means that the frequency / can be turned into a length 
L = v s /f. It might be argued that the absence of a 
quadratic term in (JjJ) prevents a straightforward defini- 
tion of such a velocity in the T — limit; nevertheless, 
it has been shown 19] that an effective phonon disper- 
sion relation at finite energy density can be evaluated 
for model J7J), yielding v s — 1.308 at e = 1. Using this 
value, we can ascertain that, at least for TV > 1000, there 
is an excellent agreement between the two approaches 
(see again Fig.OJ). 




FIG. 3: Quartic FPU model: the effective exponent a c ff of 
the finite-size conductivity for T + = 1.2, T_ = 0.8 (full dots), 
compared with the results (—S c s) of equilibrium simulations. 
The two horizontal lines correspond to the theoretical predic- 
tions, 1/3 and 2/5. 

The data presented so far rule out the value predicted 
in Ref. for the model J7J. On the other hand, such 
prediction is consistent with numerical results for hard- 
core interaction 0, . In order to test for universality 
we thus tried to bridge the two classes of models by in- 
troducing a strong repulsive potential. This has also the 
merit to remove one of the main drawbacks of a poten- 
tial like 0, namely that negative values of x n+ \ — x n 
are allowed, i.e. that particles can formally cross each 
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other if x n is interpreted as their actual position. This 
unphysical feature may somehow be circumvented by in- 
troducing explicitly a physical distance. For this purpose, 
we have added to the FPU potential a repulsive term of 
the Lennard- Jones (LJ) form 

V(z) = I(z-a) 4 + l-L-y , (10) 

where Vo is a suitable constant needed to set the mini- 
mum of the energy equal to 0. The core repulsion intro- 
duces a further time-scale, namely that of mutual "colli- 
sions" between particles |20j. Upon increasing a, at fixed 
energy, the role of the repulsive term becomes negligible 
and model l|l(J|) reduces to the purely quartic FPU J7J. 
For instance, in the region where V(z) < e the LJ energy 
contribution can be as large as 0.57 for a = 2 and e = 1, 
but it is at most 0.028, when a is increased to 2.5. Upon 
decreasing a, the LJ term progressively affects the high- 
frequency spectral range. This is because core repulsion 
becomes more relevant close to the minimum. 

In this context, one should, in principle, refer to the 
general heat-flux expression Q that, in the limit of pure 
hard-points reduces to J2 n ^n/^- Nevertheless, in the 
parameter range investigated hereby, the spectra of this 
quantity never exceed 10% of the spectra of I0J in all 
frequency channels and, more importantly, hardly show 
any singular low-frequency behaviour. We have there- 
fore kept determining the power spectrum of the flux as 
defined in Eq. (QJ. 

The effect of the LJ term on the low-frequency be- 
havior of S(f) can be appreciated already for a = 2.3. A 
direct fitting of the 3 upmost curves in Fig.^hi the avail- 
able scaling ranges) yields 5 decreasing from —0.25(0) for 
a = 2.0 up to —0.37(8) for a — 2.5. Having averaged 
the spectra over more than 5000 different samples, it is 
possible to investigate the convergence of each <5-value 
through the effective exponent JSJ). Like in the quartic 
FPU model, one can see from Fig. [21 that 5 e g decreases 
upon decreasing frequency, even on the smaller available 
range. Although on the basis of numerics alone one can- 
not exclude that the asymptotic 5- value depends on a, it 
is wiser to conjecture that the stronger the LJ term, the 
slower is the onset of the asymptotic regime. 

Altogether our simulations do not confirm the claim 
contained in Ref. @ that the hydrodynamic theory ac- 
counts for all Id models. The exponent a is found to 
be definitely larger than the expected value 1/3 and cer- 
tainly closer to the mode-coupling estimate 2/5. Any- 
how, the systematic deviations shown in Fig. [5] make 
also the convergence to this latter value somehow ques- 
tionable. In addition, we have also shown that the low- 
frequency power-law behavior is strongly influenced by 
the presence of a hard-core repulsion term: even small 
variations of the spatial scale associated with the equi- 
librium distance between interacting oscillators enhance 



finite size effects and slow down convergence with respect 
to the purely anharmonic model. This scenario rather 
suggests that at least two different universality classes 
may exist, although their physical origin is up to know 
unclear. 

A similar analysis should now be applied to the other 
models recently considered, in order to determine how 
much of the observed mutual fluctuations are the result 
of finite-size corrections. However, since we have basi- 
cally reached the limit of our computing facilities letting 
a cluster of 48 PCs run for two months, it is also clear 
that more refined analytic estimates have to be worked 
out to shed light on this puzzling scenario. 
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